Fourier Transform of the Stretched Exponential Function: Analytic Error Estimates, 
Double Exponential Transform, and Open-Source Implementation libkww. 
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An algorithm is described for computing the Laplace transform (one-sided Fourier sine and co- 
sine transform) of the stretched (or compressed) exponential function exp(— t 13 ) (also known as 
Kohlrausch-Williams- Watts function, as characteristic function of a Levy stable distribution, or 
as complementary cumulative Weibull distribution) for exponents f3 between 0.1 and 2. For low 
and high frequencies, the well-known series expansions are used; for intermediate frequencies, 
the explicit integration is strongly accelerated by the Ooura-Mori double exponential transfor- 
mation. The algorithm is implemented in C as library libkww. The source code is available at 
http : //www.messen-und-deuten. de/kww . 



PACS numbers: 07.05Kf,02.30Gp,02.30Uu,61.20.Lc 



I. INTRODUCTION 

A. Purpose 

Dynamic phenomena can be studied as function of time 
or frequency. Time-dependent relaxation functions and 
frequency-dependent spectra are related to each other 
by Fourier transforms. These transforms are only in spe- 
cial cases given by simple analytic expressions, like the 
Cauchy spectrum ( "Lorentzian" ) for exponential relax- 
ation. In general, the transforms must be computed nu- 
merically. Not rarely, the actual or perceived difficulty of 
this computation is a hurdle that prevents adequate data 
analysis. 

Specifically, relaxation in disordered systems is of- 
ten described by the stretched exponential function 
exp(— (i/r)' 5 ). However, when the same dynamics is mea- 
sured in the frequency domain, experimentalists often re- 
sort to other fit functions (like a sum of two Lorentzians 
or a Havrialiak-Negami function) that should be shaved 
away by Occam's razor since they introduce additional 
parameters without providing better fits, let alone deeper 
insights. 

This paper describes the mathematical foundations of 
a new software that allows for speedy computation of the 
Fourier transform of the stretched exponential function. 
For low and high frequencies, well-known series expan- 
sions are used. For intermediate frequencies, the Fourier 
integral is calculated explicitely. This is the bottleneck 
that determines computing times. In the present imple- 
mentation, the explicit integration is strongly accelerated 
thanks to a recent mathematical innovation, the double 
exponential transformation of Mori et al. 0, Q • 
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B. Uses of the stretched exponential 

The stretched exponential function arises in different 
mathematical contexts, for instance as Levy symmet- 
ric alpha-stable distribution, or as the complement of 
the cumulative Weibull distribution. In recent years, it 
has been found to provide a good fit to various socio- 
economic statistics, like urban agglomeration sizes, cur- 
rency exchange rate variations, or the 'success' of scien- 
tists, musicians, and Hollywood blockbusters [H, 0, HI- 

In physics, the most important use of the stretched 
exponential function is the approximative description of 
relaxation in glass-forming liquids and amorphous poly- 
mers. In 1993, Bohmer et al. Q listed stretching ex- 
ponents for over 70 materials, obtained by viscoelastic, 
calorimetric, dielectric, optical, and other linear response 
measurements. As of 2008, that review has been cited 
over 900 times, indicating a huge increase in the use of 
the stretched exponential function for describing relax- 
ation phenomena. An impressive body of literature on 
stretched exponential relaxation in molecular and elec- 
tronic glasses has been reviewed in 1996 by Phillips p}; I 
disagree, however, with his theoretical views, starting on 
the second page where he claims that j3 tends towards 1 
for high temperatures, which is not confirmed by recent 
measurements in normal and slightly supercooled liquids 
(e.g. 1,0). 

Other physical applications of the stretched exponen- 
tial function are the time dependence of luminiscence or 
fluorescence decays [l(J, and the concentration depen- 
dence of diffusion coefficients and viscosities [ll[ . In most 
applications, the exponent is restricted to values /3 < 1. 
However, in recent years some uses of the "compressed" 
or "squeezed" exponential function with 1 < < 2 have 
been proposed, mostly in protein kinetics [l2|, [l3|, [l4| , but 
also in magnetism |15| . 

The use of the Fourier transform to describe dynamic 
susceptibilities and scattering experiments has its foun- 
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dations in linear response theory. The relations between 
response functions, relaxation functions, susceptibilities, 
correlation functions, and scattering laws are briefly sum- 
marized in Appendix [A"l 



C. Historic Notes 

In physics, the earliest known use of the stretched ex- 
ponential function is by Rudolf Kohlrausch in 1854 to de- 
scribe charge relaxation in a Leiden jar. He was followed 
by his son Friedrich Kohlrausch in 1863 who used the 
stretched exponential to describe torsional relaxation in 
glass fibers, thereby improving previous studies by Wil- 
helm Weber (1841) and Rudolf Kohlrausch (1847). In the 
modern literature, these early accomplishments are often 
confounded, and a majority of references to Poggendorff's 
Annalen der Physik und Chemie is incorrect [16j |. 

In 1970, Williams and Watts introduced the Fourier 
sine and cosine transform of the stretched exponential 
function to describe dielectric response as function of fre- 
quency (l7| . Their intuition is remarkable, since they 
were neither aware of earlier uses of the stretched ex- 
ponential in the time domain, nor had they the techni- 
cal means of actually computing the Fourier transform: 
based on analytic expressions for (3 = 1 and (3 = 0.5, they 
courageously extrapolated to (3 = 0.38. In a subsequent 
paper [l8| it is quite obvious that the curves, perhaps 
drawn with a curving tool, are not really fits to the data. 

It was noticed soon that series expansions can be used 
to calculate the Fourier transform of the stretched expo- 
nential in the limit of low or high frequencies [li| [2(J Hl[ . 
Based on this, computer routines were implemented that 
complemented these series expansions by explicit inte- 
gration for intermediate frequencies [22l. |23||. In actual 
fit routines, it was found more convenient to interpolate 
between tabulated values than to calculate the Fourier 
transform explicitely [24|. Other experimentalists fit 
their data with the Havriliak-Negami function, and use 
some approximations [25l [26| to express their results as 
Kohlrausch- Williams- Watts parameters. 



D. Claims 

The present implementation improves upon previous 
work [23, [H| in the following respects: (1) The regions 
where series expansions are applicable are redetermined, 
and error estimates are corrected. (2) A wider (3 range is 
covered. (3) The integration at intermediate frequency, 
the bottleneck of previous implementations, is improved 
in speed and accuracy by using the Oouri-Mori double 
exponential transform. (4) The implementation has a 
well defined accuracy of 10~ 7 . (5) The implementation 
is made available as a C library libkww that is easy to 
retrieve, to install, and to use; the source code is released 
under the GNU General Public License so that long-term 
availability is ensured. 



Claim (2-4) require some explanation. Dishon et al. 
published tables for values of (3 between 0.01 and 2 [22| . 
However, for (3 < 0.3, these tables are useless because 
they do not cover the u> range where the Fourier trans- 
form is nontrivial. This will become clear in Sect. IrTDl 

With respect to speed of calculation, one might oppose 
that given todays computing power this is no longer a 
serious concern. However, the Fourier transform of the 
stretched exponential is often embedded in a convolution 
of a theoretical model with an instrumental resolution 
function, which in turn is embedded in a nonlinear curve 
fitting routine. In such a situation, accelerating the in- 
nermost loop is still advantageous. 

With respect to accuracy, one might argue that a nu- 
meric precision of 10 -3 or 10 -4 is largely sufficient for 
fitting spectroscopic data. However, violations of mono- 
tonicity at a level 5 can trap a fit algorithm in a haphaz- 
ard local minimum unless the minimum search step of 
the algorithm is correspondingly set to 0(5). By guar- 
anteeing full single-precision accuracy 5 = 1.2 • 10~ 7 , the 
integration of libkww into existing fit routines is made 
easier and safer. 



II. MATHEMATICAL FOUNDATIONS 

A. Notation 

We write the stretched exponential function in dimen- 
sionless form as 

f p (t) ^expH' 3 ). (1) 

Motivated by the relations between relaxation, linear re- 
sponse, and dynamic susceptibility ( Appendix |A|) . we de- 
fine the one-sided Fourier transform of fp as 

poo 

Fp(oj) := / dte iut ff,(t). (2) 
Jo 

Strictly speaking, this is a Laplace transform. It can 
be read as two-sided Fourier transform if the relaxation 
function fp is multiplied by a Heavyside step function. 
The frequency w can be construed as having a small pos- 
itive imaginary part, lo + i0 + . 

In most applications, one is interested in either the 
cosine or the sine transform, 

Q {u>) := ReF>H, 

(3) 

Vp{u) := ImFp(Lu). 

In physical applications, the stretched exponential 
function is usually written as 

f f3 At):=ex P (-(t/rf). (4) 

Its Fourier transform Fp^ T can be expressed quite simply 
by the dimensionless ([2]): 

F PtT (u) := tFp(tlu). (5) 
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B. Small u> Expansion 



with amplitudes 



For small and for large values of oj, Fp(oj) can be de- 
termined from series expansions 0, [Ifj, l2ll . HI, [23| • For 
small values of u>, one may expand the exp(iujt) term in 
([2]) . Substituting x = t 13 , and using the defining equation 
of the gamma function, 



(6) 



one obtains a Taylor series in powers of u> (in Ref. [2l( 
traced back to Cauchy 1853): 



Vp(u) 
with amplitudes 
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For > 1, these series converge for all values of w. In 
practice, however, the use of ([7]) must be restricted to 
small \u>\; for large the alternating terms become pro- 
hibitively large. 

For (3 < 1. the series ([7]) constitute asymptotic ex- 
pansions [U Hfl . In Appendix [B] it is shown that the 
modulus of the remainder is not larger than that of the 
first neglected term. This improves upon a weaker and 
unproven estimate in Ref. (23|. To minimize the trunca- 
tion error, the summation must be terminated just before 
the smallest term. 



C. Large uj Expansion 

A complementary series expansion for large u can be 
derived by expanding the exp(— t") term in @. Using 



(9) 



(where ± is given by sign a) one obtains a series in powers 
of uj-P (in Ref. [Hf attributed to Wintner 1941 [29[): 



\fc ±ikRn/2 r-, I i-l-fe/3 



Fp{w) = *£(-l) fc e 

oo 

Q/?(w) = =F (~ 1 )' C sin(fc/37r/2)S fc |w 
fe=i 

oo 



(10) 



fife = 



r(fc/? + i) 
r(fe + i) ' 



(ii) 



The convergence of these series is kind of symmetric to 
the one found in Sect. Ill Bl For j3 < 1, the series converge 
for all for (3 > 1, they (are asymptotic expansions. 

For small the are useless because the oscillating terms 
grow too much before they finally start to converge; for 
large \lu\ they converge rapidly. In the following, we will 
explicitely consider only the case oj > 0. 

To avoid numeric inaccuracies in evaluating the com- 
plex exponential or trigonometric factors for ft — > 2, 
these terms may be using the complementary exponent 



k=l 



sin(fc/37r/2)B fc w 



-l-k/3 



(12) 



and similarly for Fjg and Vp. At this occasion we note 
that there is a qualitative difference between Qp and Vp 
in the special case (3 — 2 where the large uj expansion is 
not applicable for Qp because all terms are zero. 

Concerning the computation of (jlOD for (3 > 1, two er- 
roneous staments are made in Ref. (231 ] . namely: (i) the 
most accurate results are obtained truncating the sum- 
mation before the smallest term; and (ii) the truncation 
error is less than twice the first neglected term. 

To see that statement (i) is wrong, let (3 be a "round" 
rational number like 3/2, and choose u> small so that Bk 
decreases at least for the first few values of k. Then, for 
one of these fc, the trigonometric factor in Qp(uj) or Vp(oj) 
is zero. If statement (i) were true, (JTUJ) could be truncated 
with zero error to a finite sum. This is obviously not 
true. A correct truncation criterion must be based not on 
full terms, but only on amplitudes £?&; it must disregard 
oscillating prefactors like sm(k(3w/2). 

To make sense of statement (ii), it must be reinter- 
preted as refering not to the first neglected term in its 
entirety, but only to its amplitude B n _i. But even then, 
the statement is unfounded. In Ref. [23[ it is underlaid 
by a reference to a specific page in a book on numerical 
analysis [30{ . However, that page only says "the error 
committed is usually less than twice the first neglected 
term" (the emphasis is mine), followed by a reference to a 
specific page in a 1907 book on Celestial Mechanics [3l| . 
Going back to that source, we find a rigorous theorem 
that holds only under very restrictive conditions that are 
not fulfilled here. 

A valid estimate of the truncation error has been de- 
rived by Wintner [2^| for the case (3 < 1. In Appendix [Cl 
a simplified version of his argument is presented, and 
generalized to the case (3 > 1. The result is the follow- 
ing: When the sums in (fTU)) are truncated after a term 
k = n — 1, then the modulus of the truncation error is 
not larger than 



fe=0 



(sin i 



-l-np 



B„uj 



-nP+l 



(13) 
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with 



tt/2 if/3<l, 
it/ (20) i£p>l. 



(14) 



For /3 < 1, sin0 = 1. Therefore, the truncation error is 
not larger than the first neglected term, disregarding its 
oscillating prefactor. For 1 < j3 < 2, the error bound 
increases with increasing f3. In the worst case (3 = 2, one 
has a prefactor (sin = 2 1 / 2+ ". The truncation 

criterion is: stop summation with k = n — 1 if 



(sin < 



-B„+icj 13 > B n . 



(15) 



D. Cross-over frequencies 
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FIG. 1: Cross-over frequencies ujq, 
cording to Eq. (flo) . 
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as function of f3, ac- 



The leading-order terms in ([7|) and (|10|) are power-laws 
in to. If we plot In Qp or In Vp versus In uj, these power-law 
asymptotes are straight lines. The intersection of a low-w 
and a high-w asymptote defines a cross-over frequency: 



f/3T(l + /3) sin(/?V2)\ 1/(1+/3) 







r(i//3) 

1/2 



(16) 



T(2/P) 



Fig. [T] shows these cross-over frequencies as function of j3. 
For (3 < 0.5, the curves are very steep; for (3 — > 0, both 
cross-over frequencies go rapidly to zero, the leading or- 
der being 



1//3 



(17) 



Because of this divergence, the limiting case of very small 
(3 values has no practical importance, at least in physical 
applications. 

The divergence of the cross-over frequencies (fTT)) also 
explains why previously published tables [1^ of Qp(io) 
and Vp(ijj) are useless for small exponents j3 < 0.3: As 



these tables employ a fixed, linear to grid, for /3 — > 
the nontrivial cross-over region is no longer covered; the 
tables actually contain no more than the asymptotic large 
lu power law. 



E. Numeric Integration 

Popular approaches to calculate numeric Fourier trans- 
forms include straightforward fast Fourier transform, and 
Tuck's simple "Filon-trapezoidal" rule [321] . Both meth- 
ods evaluate fpit) on an equidistant grid = kAt. The 
Filon rule optimizes the weight of the grid points. For 
small (3, the decay of fp extends over several decades. In 
that case one may either resort to brute force, increas- 
ing the number of grid points to the limits of the com- 
puters storage capacity, or use a decimation algorithm. 
The double-exponential transformation provides a sim- 
pler and more efficient alternative. Instead of optimizing 
weights, one optimizes the grid points t^, taking advan- 
tage of the zeros of the Fourier integrand. 

The double-exponential transformation was first pro- 
posed by Takahasi and Mori in 1974 for the efficient 
evaluation of integrals with end-point singularities [l|. 
Afterwards, it was improved for the evaluation of oscilla- 
tory functions by Ooura and Mori 0] • Their algorithm is 
particularly efficient for calculating Fourier transforms at 
intermediate frequencies where series expansions become 
slow or fail altogether. 

The double-exponential transform is a substitution of 
the integration variable t by a new variable k, using a 
monotonous function 4>(k) that satisfies (for the case u> > 
0) 

cj>(k -> -oo) 0, (18) 
<fi'(k — > — oo) — > double exponentially, (19) 
4>(k — > +oo) — > k double exponentially. (20) 



With 



t = -<j>{k), 



(21) 



the Fourier integral © becomes 



Fp(u) = - 



dk 



^(*)exp(i7r0(fc))/(^(fe)) . (22) 



Applying the trapezoidal rule with stepwidth 1, we ob- 
tain an approximation 

N + 

Fp(uj;N-,N+) = - ]T & fe /P) (23) 

LU — \ CO ' 

k=N- 



with scale factors 



a k = Tr<f)(k), 



(24) 
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amplitudes 

b k = 4>'{k)c k , 
and complex phase factors 

c fe = exp(iafc). 



(25) 



(26) 



These coefficients can be pre-computed independently of 
/ and u>. 

For large negative k, the conditions (fl8|) and (fT9|) en- 
sure that Ck goes rapidly towards 1 and bk towards 0. 
For large positive k, condition (|20[) implies that the ar- 
gument of the complex exponential function in (|26|) tends 
towards ink. At this point it is convenient to choose the 
summation limits according to 

{half-integer when calculating Qa, 
P (27) 
integer when calculating Vp. 

In consequence, when calculating the cosine transform 
Qp, all k are half-integers, and therefore Re Cfc goes to 
0. Similarly, for the sine transform Vp, all k are integers, 
and Im Ck — > 0. So, in any case, the relevant part of Ck 
tends rapidly towards as k — > ±oo. This allows us to 
approximate (j2"2"|) by (f2"3")) with surprisingly low values of 
|JV±|. 

All transformations considered by Ooura and Mori [2| 
have the form 



(28) 



1 - exp (-x(fc)) ' 
The function x(fc) must fulfill the conditions 

x(k — ► — oo) — ► — oo exponentially, (29) 
X (0) = 0, (30) 
x(k —* +oo) — > oo exponentially, (31) 



Condition (|30|) guarantees that numerator and denomi- 
nator of l|28p have a zero at the same location k = so 
that the singularity is removable, 



The slope at this point is 
*'(0) = ; 



1 



X'(0) 



x"(Q) 

2x'(C)) 2 



(32) 



(33) 



These values are needed to compute the coefficients ao, 
bo for use in the sine transform. 

We note in passing that (|2"6")) becomes problematic for 
large positive k. Following Ref. Q, we use ([28]) to rewrite 
c k as 



Ck = exp(i7rfc) exp 



iirk 



e x(fc) - l 



(34) 



For the calculation of Q/3, one has to take the real part of 
(fM)) . with k being half-integer values; for Vp, the imagi- 
nary part, with integer k. The resulting expressions can 
both be written in the same form, using the floor function 
Lfcj: 



Re Ck 
Im Ck 



= (-1 



JfeJ 



irk 



e x(fc) 



(35) 



To proceed, the function x(fc) must be specified. Orig- 
inally, Ooura and Mori [33| had proposed 



Xi(k) = 6 sinh(/ifc) 



(36) 



where h is a parameter that controls the accuracy of the 
trapezoidal approximation. Later, after studying the in- 
fluence of singularities near the real axis, they suggested 
1 



Xa(*) = 8M; 



47h 



with 



lh = v/l + ln(l + 7r//i)/(4/i). 



(37) 



(38) 



Since the stretched exponential function has no such sin- 
gularities, there is no reason to prefer (|3"T)) over (|3^|) . 
Numeric experimentation shows that the integration is 
rather insensitive to the choice of x- 



III. IMPLEMENTATION 

A. Accuracy 

Let as write S for either Qp(uj) or Vp(uj). The numeric 
computation of S shall aim at a relative accuracy 5. We 
approximate S by the sums (O , (|10p , and (f23|) , which we 
denote for short as 



5n 



n-l 

E 



Sfe. 



(39) 



In the cases of the small and large ui expansions ([7|) 
and (HHJ), we have derived in the appendices [B] and [Cl 
upper bounds for the truncation error: 



\S n — S\ < e n . 



(40) 



The sum (1391) is computed incrementally until either the 
relative error e n /S n is smaller than 5 or the computation 
fails for one of the reasons to be discussed in a moment. 

The accuracy of the trapezoidal approximation (|23p 
depends on the step width, expressed by the param- 
eter h in (|36] ) or ((37 |) . For given ft,, N_ and N + 
(and thereby n = N + — N_ + 1) are chosen such that 
\bj^ ± \ < 5. The trapezoidal sum SVi,7i is calculated for 
decreasing step widths hj until either the relative change 
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\{Sn,hj - Sn^^J/Sn^] is smaller than S or the compu- 
tation fails. 

Computations can fail at any k < n for one of the 
following reasons: 

(i) Sfc has grown beyond the largest floating-point num- 
ber; 

(ii) Sk has decreased below the smallest normalized 
floating-point number; 

(iii) in the case of an asymptotic series, eu+i > e^: the 
terms start to grow before the desired acccuracy e n /S n < 
S has been reached. 

Furthermore, computations can be declared ex post to 
have failed: 

(iv) While incrementing ((39)) , we also update a variable 
that contains the largest amplitude contributing to S n , 
z n = maxfc< n |sfe|. If z n ^> S n , then the accuracy of S n 
suffers from the cancellation of huge alternating terms. 
Let e denote the accuracy of the floating-point data type 
used for computing the Sfc and S n . Then the accuracy of 
S n is at best ez n . We only accept S n if ez n < SS n . 

In our implementation, we aim at single-precision ac- 
curacy 5 = 2~ 23 ~ 1.2 • 10~ 7 . Internally, Sk and S n are 
double-precision variables with e = 2 -52 ~ 2.2 • 10~ 16 . 



B. Results and Ranges 
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FIG. 3: Numeric Fourier transforms as in Fig. [2] but for a 
much narrower frequency range. The dashed line for (3 — 2 is 
the Gaussian obtained by analytic cosine transform. 



/?, both expansions fail for one of the four reasons enu- 
merated in Sect. IIII AI In these cases, one must resort 
to numeric quadrature based on the double exponential 
transform (solid lines). 
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FIG. 2: Numeric Fourier transforms Qp(iv), Vp(ui) over wide 
logarithmic ranges. Discrete symbols are obtained from the 
expansions (O for small u> (blue) and (|10l) for large ui (red); 
solid lines are results of numeric quadrature using the double 
exponential transform. 

Fig. [2] gives a broad overview over numeric results; 
Fig. [3] concentrates on the crossover region between the 
two asymptotic power laws. Blue symbols have been ob- 
tained by the small u> expansion ([7]) , red symbols by the 
large lo expansion (|10p . For some combinations of u> and 



0.8 1.2 1.6 

P 

FIG. 4: Limits for the two cosine transform expansions l[7| 
(blue) and (|10[) (red) in the j3, uj plane. The lower frame 
contains an enlarged view for /3 ~ 1. 

Fig. 0] shows the large w limit of the small to expansion 
and the small uj limit of the large u> expansion for the 
case of the cosine transform as function of j3. Results for 
the sine transform are very similar. These limits have 
been calculated by a simple interval division alogrithm, 
using a Ruby script kwwlimits . rb that is part of the 
source distribution. For intermediate values of /3, there 
is only a tiny gap between the two expansions; for the 
sine transform near f3 ~ 1.4, they even overlap. On the 
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other hand, for (3 < 0.3 the gap between the two series 
expansions extends over more than a decade in to. 

The data points used in Fig. |4] are copied from the 
output of kwwlimits . rb into the source of the high-level 
routines that implement Qp(uj) and Vp{io) for all values 
of lu and for 0.1 < (3 < 2. For given (3, m a simple inter- 
polation is used to decide whether a series expansion is 
likely to work or not. In rare cases (for to very close to a 
border line) , the decision will be wrong. In case of a false 
negative, some computation time is wasted doing the nu- 
meric quadrature though the series expansion would have 
worked. In case of a false positive, the subroutine that 
implements the series expansion will return an error code; 
then the high-level routine will resort to calling the sub- 
routine that implements the numeric quadrature. 

C. Special Case /3 —> 2 

The cosine transform requires special care in the case 
j3 — * 2. In the limit (3 = 2, the transform is just a 
Gaussian, 

Q 2 H = ^exp(-c 2 /4). (41) 

For (3 < 2, Qp is very close to Q2 for low frequencies, 
whereas for high frequencies it has a power-law tail 

Q 2 _p{uj) ~ sin(/37r/2)r(2 - (3)uj- 3+ ^ for to > 1 (42) 

that vanishes for (3 — * 0. 

In the cross-over range between these two regimes, the 
numeric quadrature fails to reach the required accuracy 
because of cancellation. This problem can be solved by 
quadrating not fp(t) but the difference fpit) — fz(t). The 
analytic transform Q2(w) is then simply added to the 
numeric Re FT[f@ — / 2 j. 

D. Public Interface 

Routines for the computation of Qp(uj) and Vp{u)) 
have been implemented in form of a small library 
libkww. In order to ensure maximum portability, the 
programming language C has been chosen. The source 
code is published under the terms of the GNU Gen- 
eral Public License (GPL) [34]. The project website 
is http://www.messen-und-deuten.de/kww. The entry 
page summarizes the scope of libkww in the style of a 
Unix manual page. A link leads to the source code direc- 
tory. 

To prepare the source code for distribution, standard 
procedures of the free software project GNU have been 
used. After downloading a compressed source package, 
the standard command sequence 

$ tar zxvf kww- <v ers i on> .tgz 
$ cd kww-<version> 



$ ./configure 
$ make 

$ sudo make install 

is sufficient to install libkww on a Unix system. 

Following Unix manual page conventions, the applica- 
tion programming interface (API) can be documented by 
the "synopsis" 

ttinclude <kww.h> 

float kwwcf (float omega, float beta); 
float kwwsf (float omega, float beta); 

The letters c and s stand for cosine and sine transform, 
respectively; the routines return Qp{w) and Vp{uj). Fol- 
lowing a convention of the C standard library libc, the 
letter f stands for the data type float of the return 
value. By choosing this single-precision data type instead 
of the more common double-precision type double, it is 
made clear to the user that no attempt has been made 
to achieve more than single-precision accuracy. 
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APPENDIX A: DESCRIPTION OF RELAXATION 
IN TIME AND FREQUENCY 

The use of the Fourier transform to describe dynamic 
susceptibilities and scattering experiments has its foun- 
dations in linear response theory. In this appendix, the 
relations between response functions, relaxation func- 
tions, susceptibilities, correlation functions, and scatter- 
ing laws shall be briefly summarized. 

The linear response B(t) to a perturbation A(t) can be 
written as 

B(t)= [ dt' R(t-t')A(t'). (Al) 
J —00 

Consider first the momentary perturbation A(t) = 6(t). 
The response is B(t) = R(t). Therefore, the memory 
kernel R is identified as the response function. 

Consider next a perturbation A(t) = e vt Q(— f) that is 
slowly switched on and suddenly switched off (O is the 
Heavyside step function, r\ is sent to + at the end of 
the calculation). For t > 0, one obtains B(t) = <I>(i) 
where $ is the negative primitive of the response func- 
tion, R(t) = —dt^if). Since $ describes the time evo- 
lution after an external perturbation has been switched 
off, it is called the relaxation function. Kohlrausch's 
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stretched exponential function is a frequently used ap- 
proximation for $(£). 

Consider finally a periodic perturbation that is 
switched on adiabatically, A(t) = exp(—iu->i + rjt), imply- 
ing again the limit 77 — > + . The response can be written 
B{t) = x(w)A(t), introducing a dynamic susceptibility 



dte l(uJ+lv)t R(t). 



(A2) 



This motivates our definition Q of the one-sided Fourier 
transform F(u) of the relaxation function $(f). Since 
R(t) = —dt$>{t), there is a simple relation between x and 
F: 



X(v) = *(0) + iuF p (u). 



(A3) 



In consequence, the imaginary part of the susceptibility, 
which typically describes the loss peak in a spectroscopic 
experiment, is given by the real part of the Fourier trans- 
form of the relaxation function, Im x = wRe Fp(uj). 

Up to this point, the only physical input has been 
Eq. (|Alj) . To make a connection with correlation func- 
tions, more substantial input is needed. Using the full 
apparatus of statistical mechanics (Poisson brackets, Li- 
ouville equation, Boltzmann distribution, Yvon's theo- 
rem), it is found |35[ that for classical systems 



(A(t)B(0)) = k B T$(t). 



(A4) 



Pair correlation functions are typically measured in scat- 
tering experiments. For instance, inelastic neutron scat- 
tering at wavenumber q measures the scattering law 
S(q,ui), which is the Fourier transform of the density 
correlation function, 



S{q,u) 



1 

2^ 



die*"*<p(g,t)V(g,0)>. 



(A5) 



In contrast to (JSJ) and (|A2|) . this is a normal, two- 
sided Fourier transform. In consequence, if we let 
(p(q,t)*p(q, 0)) = <&(£), then the scattering law S(q,u) 
is proportional to the real part Re F(cu) of the one-sided 
Fourier transform of <&(£). 



APPENDIX B: TRUNCATION ERROR IN 
SMALL u EXPANSION 

In this appendix, we derive an upper bound for the 
error made by truncating the small u> expansion (J7|| . We 
consider the cosine transform, and write the Taylor ex- 
pansion with Lagrange remainder as 



n—X k 

fc! 

k=0 



Q^M = IJQ^ (0)^- + 0^(0^- (Bi) 



ni 



with < £ < w. From 0, we know that 

f for fc odd, 

I (—1) Afc for fc even. 



(B2) 



Choosing n even, we have 
\n( n h 



= 1 


Re F< n) 


(01 


< 1 








-£r- / die^e" 4 " 




/>OC 

/ di(it)"e 4 «*e- t ' 3 
Jo 


J 


poo 

/ dt 




iitfe^e- 1 " 


J 


/>00 

/ dit 
'0 






ff>(0) 






Oi n) (o) 





(B3) 



Therefore, the truncation error is not larger than the first 
neglected term. The same holds for the sine transform. 



APPENDIX C: TRUNCATION ERROR IN 
LARGE u EXPANSION 

In this appendix, we derive an upper bound for the er- 
ror made by truncating the high-w expansion (|10p . sim- 
plifying an argument of Wintner [2^ | , and generalizing it 
to cover not only the convergent case f3 < 1 but also the 
asymptotic expansion for /3 > 1. Substituting s = cut, 
the Fourier integral Q takes the form 



P 00 

Fp(uo) = uj^ 1 I ds exp (is — oj _/3 s' 3 ) . 
Jo 



(CI) 



For brevity, we discuss only the cosine transform 
Qp(w) = Re Fp(u}), which we rewrite as Qp(u) = 
G{uj~P)/lu, introducing the functions 

/>oo 

G{x) = Re / dsj(s,x,0) (C2) 
Jo 



and 



7(s, x, a) — s a exp (is — xs 13 ) . (C3) 



The Taylor expansion of G(x), including the Lagrange 
remainder, reads 



G( a; )=^GW(0) ¥ + GW(0- [ (C4) 

fc=0 

with < £ < s and 

/>oo 

G« (£) = (-!)* Re / ds 7 (s,£,fc/3). (C5) 
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Now, we choose an integration path C in the complex 
plane, consisting of two line segments, s and se 1 ^ , and 
two arcs, re lip and Re l(p , with 0<r<s<R<oo and 
< <p < <j>. The integral of 7 along this path is zero: 



dz 7(2;, x, a) = 0. 



(C6) 



The contributions of the two arcs tend to as r — > 
and R — > 00. Hence the contributions of the two line 
segments have equal modulus. This allows us to obtain 
the following bounds: 



\G (n) (0\ = 



< 



< 



(-l)"Re / ds 7 (s,$,n[3) 




ds7(s,^,n/?) 



00 

ds7(se^,£,n/3) 



00 

ds |7(se^,^,n/3) 



/ ds I s «/y* n/3 exp (ise^ - ^e 

/■oc 

/ ds s n/3 exp (— ssin(£ - cos(</>/3)) 
Jo 



(C7) 



At this point we choose 

fvr/2 if/3<l, 

= 

U/(2/?) if/3>l, 
which ensures cos(</>/3) > 0, yielding a bound 



(C8) 



G {n) (0 



< / dss n/3 



s ^ exp s sin < 



(C9) 



that is independent of £. Evaluating the well-known in- 
tegral ([H]) we obtain 



< 



r(n/? + i) 

(sin 0) 71/3+1 ' 



(CIO) 



Only trivial changes are needed to adapt this argument 
to the sin trafo Vp(u>). 
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